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Abstract 

In this paper, we study the numerical simulations for Euler system 
with maximal density constraint. This model is developed in |9l [TT] 
with the constraint introduced into the system by a singular pressure 
law, which causes the transition of different asymptotic dynamics between 
different regions. To overcome these difficulties, we adapt and implement 
two asymptotic preserving (AP) schemes originally designed for low Mach 
number limit [161 118] to our model. These schemes work for the differ- 
ent dynamics and capture the transitions well. Several numerical tests 
both in one dimensional and two dimensional cases are carried out for our 
schemes. 

Key words: Finite volume scheme. Congestion, Asymptotic-Preserving schemes, 
All-speed flows, Pressureless Gas Dynamics 

1 Introduction 

Several models involve congestion constraints: concentration constraints occur 
in multi-phase flow modeling |12| , maximal density constraints occur when deal- 
ing with finite-size interactive agents in herds of gregarious mammals |17j . in 
cars or pedestrians flows [HllinilS], flux constraints occur for supply chains 
The dynamics in congested regions strongly differ from the dynamics dynamics 
in free regions. To study the transitions between congested and free regions, 
a general methodology was first carried out in [12] for multiphase flows and 
later on generalized to traffic [9, or herding problems [T7]: the stiffness of the 
constraint leads to a singular perturbation problem and then the limit prob- 
lem provides a clear cut-off between the two dynamics. In this paper, we will 
consider the Euler system with a singular pressure which encodes a maximal 
density constraint. As in |17j . the limit problem is a two-phase model between 
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incompressible regions, where the maximal density is reached, and compressible 
regions for lower densities. Our goal is to provide numerical schemes that are 
able to capture this limit problem and these transitions. In this paper, we adapt 
and compare two numerical methods presented in |18j and in |16j for the low 
Mach number limit. 

A lot of efforts have been made to devise numerical schemes valid for all 
Mach numbers, that is, for both compressible and incompressible flows. They 
avoid the switch between different methods, when different Mach numbers occur 
in different sub-domains. Among such schemes, one approach is the extension 
of compressible conservative methods to incompressible flows thanks to precon- 
ditioning techniques [371 HOI UHl HE] • The second approach is the extension of 
incompressible methods to compressible flows and leads to pressure correction 
methods on staggered grids flTj and their conservative versions [33l[6l[38]. Their 
adaptation to the conservative frameworks has led to time semi-implicit schemes 
: the implicit discretization of the fluxes (the mass flux and the pressure part 
of the momentum flux are taken implicitly) is combined with the resolution of 
the elliptic equation satisfied by the pressure. The implicit treatment of the 
pressure flux ensures stability with respect to the propagation of fast acoustic 
waves in the low-Mach number limit but induces a lot of diffusion. We can cite 
numerous works following this methodology [551 IMl 1301 1311 131| • The methods 
we consider in this paper are among the simplest ones: the scheme in [TS] is a 
semi-implicit scheme with a division of the pressure into explicit and implicit 
parts and in UJJ , the Gauge decomposition of the momentum enables the hy- 
drostatic pressure to act only on the divergence-free part of the momentum. 
The former method will be called in the present paper the Direct method, while 
the Gauge method will refer to the latter. 

The purpose of this paper is to present simple variants of the Direct [TH] 
and the Gauge method [TB], that are able to handle congestion problems. As 
announced above, they are designed to solve the isentropic Euler system sup- 
plemented by a pressure law p{p), which is singular as the density p approaches 
a maximal density denoted p* . A small parameter e is introduced to control 
the stiffness of this maximal density constraint: the rescaled pressure ep{p) is 
of order 0(f) in congested regions p ^ p* and of order 0(e) in low density 
regions p < p* . In the limit £ ^ 0, the system leads to a two-phase model: 
the incompressible Euler system in maximal density domains and the pressure- 
less gas dynamics system for uncongested densities domains. This asymptotic 
model was first proposed and studied in [121 [5] in a one-dimensional framework. 
However, this asymptotic model is only partially defined since transmission 
conditions at the interfaces between the two phases are lacking. Besides, un- 
less one-dimensional solutions can be provided (see |I2j and appendix [A|) , their 
extensions to the two-dimensional case are open problems (especially the dynam- 
ics of two colliding congested domains). In this context, asymptotic preserving 
scheme is a good tool for this problem. 

The Direct and the Gauge methods are called asymptotic-preserving (AP) 
since they are uniformly consistent with the low-Mach number limit. Besides, 
they are also uniformly stable. These methods are expected to capture both the 
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compressible and the incompressible dynamics arising in the congestion limit 
of the Euler system with the maximal density constraint. In this case, such 
AP numerical schemes are very powerful since they provide the dynamics of 
transitions, for which analytical results may be lacking. Moreover, they enable 
us to avoid dealing with physical and numerical interface tools, such as front- 
tracking [35] or volume-of-fluid methods [S^. Unlike these tracking methods, 
ours are front-capturing methods and then share some analogies with level-set 
[32] and diffusive interface methods [1]: like level-set method, the dynamics of 
the transition are implicitly embodied in the dynamics of an auxiliary function 
which here is the density and like the diffusive interface methods, the sharp 
interface is viewed as the limit of the smooth transitions of the perturbation 
problem. 

The Direct method cannot be directly applied to the congestion problem. 
Indeed, in [18], the pressure p{p) is splitted into an explicit part Pq{p) and an 
implicit part pi (p) in order to keep some numerical diffusion and avoid numerical 
oscillations. For the singular congestion pressure-law, we modify this splitting, 
such that it still ensures the stability of the scheme. Besides, it ensures the 
consistency of the explicit part of the scheme with the limit pressureless gas 
dynamics in the low density regions. Indeed, the pressureless gas dynamics 
system is weakly hyperbolic and there is no uniqueness of the entropic solution 
[TT] . Then, keeping an explicit pressure po (p) makes the asymptotic numerical 
solutions consistent with the good entropic solutions. For the same reasons, this 
pressure splitting is introduced into the Gauge method. 

The AP property of the two methods for the congested Euler system is 
demonstrated for the congested domains. This analysis is hard to extend to co- 
existent congested and uncongested regions since the dynamics of the interfaces 
between the different regions is not explicitly implemented into the schemes. 
However, several numerical test-cases provides numerical evidence of the AP 
property. Comparisons of the two schemes are also carried out and different 
behaviours of the schemes at the interfaces are measured. 

The paper is organized as follows. In section 2, we introduce the Euler 
system with the maximal density constraint. To have some basic idea of the 
solution, we give its formal asymptotic limit. In section 3, we describe the 
time semi-discretization of the Direct and Gauge schemes. The fully space and 
time discretizations are exposed in section 4, in one dimensional setting: they 
are based on the centered Rusanov scheme, also called the local Lax-Friedrichs 
scheme [27] . The two dimensional setting is quite similar. Therefore, we present 
it in the appendix. Finally, numerical simulations are performed in section 5 to 
compare the two schemes: several test cases in the one dimensional setting and 
one case in the two dimensional setting. Two appendices close this study: we 
describe one-dimensional solutions and the two dimensional discretizations. 
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2 The Euler system with congestion and its asymp- 
totic Hmit 



2.1 The model 

We consider the two-dimensional Euler system: 

dtP + ^^-q^ 0, (1) 
9tg + V,- (^^^ +V,p(p)=0, (2) 

where p(a;, e M denotes the mass density, q — pu{x, t) £ R'^ is the momentum 
field depending on the position x gM."^ and the time t > 0. The pressure p{p) is 
an increasing function such that p{p) ~ p'^ for densities p <^ 1 and p{p) — > +00 
as p tends to the congestion density p* . In the following, we will consider the 
function: 

P(P) = 7 'T>0. (3) 



p* , 

This pressure prevents the density from exceeding the congestion density. A 
variant is the van der Waals equation of state |22]. The operators Vx and V^;- 
are the gradient and the divergence of vector fields or tensor. For two vectors 
a and 6, a (g) b denotes the tensor product. 

This model already appears in [T7] with the additional constraint: q/p — u £ 
S^. In this paper, we focus on the pressure singularity and the corresponding 
numerical schemes. 

The singular pressure induces two different dynamics: for regions with densi- 
ties near p* , the pressure takes very large values in comparison with the pressure 
in low-density regions. As in [17] and previous work on traffic modeling we 
would like to clearly separate the two different dynamics. To this aim, we rescale 
p{p) into sp{p) , where the parameter e ^ 1 is the scale of the pressure in the low 
density regions: p{p) = 0{e) for density p <C 1 while p{p) = 0{l) for density 
p - p*, see Fig.[T] 

Denoting and q' the new unknowns, system ([T])-([2]) becomes: 

dtp' + V, • q' - 0, (4) 

Moreover, taking the time derivative of Q and subtracting the divergence of 
equation ([s]), we easily obtain the wave-like equation satisfied by the density: 

,2 . i'q'®q' 



dt-P - : y^-^j - = 0, 

where denotes the tensor of the second derivatives and for two tensors a 
and b and a : b denotes the contracted product of tensor. The operator 
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Figure 1 : The pressure p and the rescaled pressure ep 

denotes the Laplacian of a scalar function. Actually, system (|4])-([5]) is a strictly 
hyperbolic problem, with characteristic wave speeds in the a;-direction (where 
X is the first component of a basis of M^) given by: 



(6) 



where w| is the a;-component of the macroscopic velocity u'^ (x, t) — q'^{x, t) / p'^ (x, t) . 
Standard hyperbolic numerical schemes have to resolve the Courant-Friedrichs- 
Levy (CFL) condition: 



max(|Af|,|A^UA||)Ai< Ax. 



(7) 



In the next section, we will see that this constraint may be too stringent for 
these schemes to capture the asymptotic limit. 



2.2 The asymptotic limit 

The limit of the pressure term sp{p^) depends on the limit of p"^ . Indeed, if 
p^ ^ p with p < p* , then £p{p^) converges to 0. Otherwise, p^ — p* and the 
limit of ep{p'^), denoted p, can be non zero and depends on the convergence rate 
of p^. We assume that the limit p is always finite. Therefore, the formal limit 
of system Q-ljsj) as e goes to zero is: 



dtp- 
dtq 



q-0, 

'q® q 
P 



V^p = 0, 



{p-p*)p^O. 



(8) 
(9) 
(10) 



A one-dimensional version of this asymptotic model was proposed for two-phase 
flow modeling in [1^, where the density plays the role of the volume fraction 



5 



of liquid in a liquid-gas model. The derivation of the model lies on a relaxation 
to zero of the relative velocities of the gas and liquid and is therefore different 
from the one studied in this article. Existence and stability of solutions are 
proved for the one-dimensional version of system ([8])-([9])-( 10 1 in [5], and in all 
dimensions with viscous term in [29] . 

As regards the characteristic speeds, we note that if ep{p^) tends to p < -l-oo, 
then we have p* — = O(e^) and then ep'{p'^) = 0{e^~^). Therefore, if 
7 > 1, Aj_ can become infinite and waves with infinite speed can occur. It 
is the low-Mach number asymptotics that leads to incompressible dynamics. 
Actually, in the congested domain where p — p* , system ([8|-([9|-( 10 ) yields the 
incompressible Euler equation: 

P = P , 

V^-u = 0, (11) 

dtU + U ■ Vajtt + — Va;P ~ 0. 

P* 



Equation (11) is the incompressible constraint and the Lagrange multiplier re- 
lated to this constraint is the pressure p. The CFL condition ([t]) degenerate 
into At = 0: standard hyperbolic schemes are unable to compute the asymp- 
totic dynamics. Thus, numerical schemes with relaxed CFL condition have to 
be designed. 

In the free domain where p < p* , the CFL condition Q is not an obsta- 
cle although the system degenerates into a non-hyperbolic problem: limAf = 
lim A3 = u. This is a large-Mach number asymptotic. Numerical schemes, 
originally developed for hyperbolic systems, have to be proved to capture this 
singular limit. System ([8])-([9])-( 10 1 yields the pressureless gas dynamics: 



P<P , 

dtp + ^x ■ 0, 

a., + v..(^)=o. 

Without upper-bound on the density, this system would lead to concentration 
phenomena even for smooth initial data and so the density may become a mea- 
sure with singular part. It is related to the so-called sticky particle dynamics. 
Existence of solutions and numerical schemes have been developed in the one- 
dimensional case [m [T3] . 

The asymptotic system ([8])-([9])-( 10 1 is not complete: the well-posed defini- 



tion of p requires boundary conditions at the interface between congested regions 
{p = p*} and free regions {p < p*}. They can not be directly obtained by the 
formal derivation. One possible answer to this question is to look at the theo- 
retical asymptotic behaviour of solutions of the initial system ([4])-([5]). Such an 
approach is investigated in appendix [X] but only in the onc-dimensional case. 
Extensions to two dimensional settings are difficult and will be the subject of 
future works. The second possible approach is to use numerical schemes to 
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capture the asymptotic dynamics: this is the main methodology we develop in 
this paper. It has the advantage to be applicable in any dimensions. How- 
ever, one-dimensional problems are still good test cases to valid these numerical 
schemes. 



3 Time semi-discretization schemes 

This section is the center of this paper. It is dedicated to the presentation of two 
numerical schemes, which are able to capture the asymptotic limit of the Euler 
system with congestion presented in the previous section. For this purpose, we 
adapt the asymptotic-preserving (AP) schemes developed in JB] and [TB] for 
the low-Mach limit of the isentropic Euler system. 

3.1 The time semi-implicit discretization 

We first define a time semi-implicit discretization, which will be the building 
block of the considered AP schemes. 

Let p",(7" be the approximations of the density p and the momentum q at 
t" = nAt, n — 0, 1, . . ., where At is the time step. The semi-discretization of 
the AP scheme for the n-th time step is as follows: 



At 



V, • = 0, (12) 
+ V.Mp"+1))=0. (13) 



At 

The full discretization in time and space is postponed to the next section. We 
want to show that the scheme is asymptotic-preserving. In other words, it 
captures the correct behaviour of the limiting equation as e — > 0. To achieve 
this, the implicitness of p, q is crucial. 

Observe that the explicit part of the above scheme is pressureless. However, 
the pressureless Euler system is weakly hyperbolic, giving rise to the formation 
of density concentrations known as delta-shocks. Several numerical schemes 
for this system was proposed in [T^l US IH [7] . In [T3] , the authors proposed 
a kinetic scheme, that is valid for the isothermal Euler system and leads to a 
kinetic scheme for the pressureless system in the vanishing pressure limit. Here, 
to avoid this difficulty and as already proposed in [TS] , we split the pressure into 
an explicit and an implicit part. Numerical tests in section [5] will demonstrate 
how this splitting reduces oscillations. Thus, the scheme is written as follows: 

„n+l _ n 

P P j_ Y7 „n+l 



At 



V, • = 0, (14) 

V. • + V,(epo(p")) + V.(£pi(p"+i)) = 0, (15) 
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Figure 2: The plots of ep and epo as functions of p with e — 10~^, 7 = 2. 



where the exphcit part is given as 



ri 



Poip) = < 



Mp* - 5) + p' [p, - 5){p - + 5) 
+ \p"{p*-5){p-p.+5f) 



p < p* - S, 



a p > p^ ~ s, 



(16) 



and the imphcit part is 



Pi{p) ^ pip) - Pa{p), S = e-'-* 



(17) 



The choice of 5 makes sure that po and its derivatives up to second order are 
always bounded. To make sure all the coefficients appearing in the elliptic 
equation we will derive in the next section are continuous, we choose po to 
be a second order approximation to p, instead of a first order one. For later 
usage, also note that the function pi{p) is invertible. This is easily seen from 
the property of function p and pq, see Fig[2j By the definition, the Courant- 
Friedrich-Lewy (CFL) condition for the explicit part is 



At < 



aAx 



^{\f,\ + vwp)y 



(18) 



where a is the Courant number. Since epQ is always bounded, the CFL condition 
can be satisfied uniformly in e. 
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3.2 The Direct method 

To get the solution, we will rewrite the above scheme into another form. By 
inserting (15) into (14 1, we can get an elliptic equation for p: 

^ ^/ + V, • q" - AtA,(epi(p"+i)) - Atv2 : [^-^^ ~ AtA, (epo(p")) 

(19) 

From this equation, we can solve However, if we solve p directly, the 

density constraint p < p* may not be satisfied due to discretization errors. 
Thus, we write — p{pi^^) in (19 1 and solve the equation in terms of pi. 

The density constraint p < p* will be automatically satisfied. Moreover, the 
positivity of p can be ensured by the fact that the discretized equation satisfies 
the maximal principle. 

Once p""*"^ is obtained, we can obtain q"+^ from ( 15 ) easily. 

g„+i ^ g„ _ . (^^1^^ + V, (epoipn) + V,(epi(p"+i))| . (20) 

Remark 1 In the numerical simulation, we can also improve the accuracy by 
implementing a fully implicit scheme, which iterates the above scheme to solve 
(14) and (15) with implicit. Suppose p"-+^'^ = and = q" and 

pn+i,k ^ qTi+i,fc solutions to the following equations. 

n+l,/c + l _ n 

P + • q" - AtA^(epi(p"+i''^-+i)) 

/ „n+l,k ^ „n+l,fc\ 

-AtVl: i^^^—;^, j + A. (epo(p")) = 0, (21) 

gU+l.k+l ^ g„ _ [^—Sii j + i^Poip")) + V.(£pi(p"+l'^-+l)) 

(22) 

As k — > oo, t/ie solution approximates to the one solving the fully implicit scheme 



P 

the additional computational cost 



(both in p and ^^J. This modification provides little improvement compared to 



3.3 The Gauge method 

Another way to implement the AP scheme is the Gauge method developed in 
|16j . It can be obtained by applying the Gauge decomposition 

q = a V^Lp, • a = (23) 

where a is the incompressible part of field q and (p is the irrotational one. This 
decomposition is expected to be more robust for capturing incompressibility 
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constraint. We will see that it is partially right. By including this decomposition 
into equations (14) and (15 1, we get 



At 



AtA,(£Pi(p"+i)) 
(g) q 



AtV 



AiA,(epo(p")) - 0, 



A..P"+i 



1 

At 



,n+l 



a' 



At 



A.(epo(p")), 
-V,(£Po(p")) 



V,P"+i =0, 



n+l 



(24) 

(25) 
(26) 

(27) 
(28) 



Indeed, the equation (24) for p"^^ is derived from (14) and (15) similarly as in 
the Direct method. The Laplace equation ( 25 ) for (^"^^ is the direct consequence 
of applying the decomposition ( 28 ) and Va, • a"+^ = to the density equation 
(14). The equation (26) for P"^^^^nd the equation (27) for a"+^ are obtained 



by inserting ( 28 1 and Vx • a" 



Va; • a" = into the momentum equation 



(15). Here a new unknown P is introduced, which is defined by 



p 



.ri+1 _ 



if" 



At 



(29) 



since (I24|)-(l25j) and ^ imply (jlSj). 

The original equations (14) and (15) can also be recovered from (24)-(28l 
by assuming V^, • a" = 0. In fact, the equations for P"+i and a"+^ ((26)-(27)) 
and Vx • a" = imply V^; • a"+^ — 0. This leads to the density equation (14 1 
from the i^""*"^ equation (25) and the decomposition (28). (14) combined with 
the Pi^^ equation (24) will then allow us to recover the momentum equation 



The boundary conditions By solving the equations (24)-(28) in sequential 
order, we can update the value of p, q. To do this, we need to provide boundary 
conditions for the Laplace equations for Lp and P. The boundary conditions for 



pn+i g^j^g somehow straightforward due to the implicit relation (29), once Lp is 
known. And in solving (p""^^ , we may impose Dirichlet boundary condition on 
(/j""*"-^ as follows: 

^"+'|n = 0. (30) 

Indeed, other non-homogeneous Dirichlet boundary conditions can be chosen. 
This makes the unknown determined up to a linear function in space. 

However, this uncertainty can be removed by redefining 0"+^. So we can always 
impose the homogeneous Dirichlet boundary conditions. 
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Simplification in the one dimensional case Observe that in the one di- 
mensional case, V X -0 = implies that a is independent of x. We may thus 



rewrite ( 26 ) and ( |27[ ) into a simpler equation by using ( 29 1 



At 



+ £Po(p")+epi(p"+') 



1 



(31) 



where the space-domain is [b, c] and f\l — f{c) — f{b). This equation can be fur- 
ther simplified, since we impose the homogeneous Dirichlet boundary condition 
on tp all the time. This lead to 



a"+i = a" - 



At ({q^) 



P" 



+ £po(p") + £Pi(p"+') 



(32) 



With this reformulation, in the one dimensional case, we can reduce the number 
of elliptic equations to be solved to one and update the space independent 
variable a more efficiently. Also as a consequence of (30), a" should be deffiied 
as the average of g". 



1 



q^dx. 



(33) 



In summary, the Gauge method in the one dimensional case is implemented 
through equations pll), psl), 1281) and (l32l). 



3.4 Discussion of the AP property 

As for the AP property of the scheme, we give a formal proof. To be precise, we 



want to show that the system (19 1 and (15) becomes the incompressible Euler 
system as e — in the congested region. 

However, in contrast with the low Mach number limit of the isentropic Euler 
equation discussed in [18^ or [16J, the singularity in our model is embedded in the 
definition of p. New congestion regions may arise from non-congested ones. And 
by the analysis in section]^ there is the possibility that p'^ p* but the limit 
of the pressure epi(p^) — ^ as e — > 0. This means that although the congestion 
density is reached, there is no real congestion in this region. So it seems better 
to characterize the congestion region by defining p = lim£_j.o ^Pi{p^) > 0. 

Currently, we can only show that in regions where both p"^^ = lim£_j.o epi ((p"^^)^) > 



and — lime_j.o epi((p")'^) > 0, (19) and (pO| tend to the incompressible 



Euler system as e — 0. The assumption lim^^Q epi{{p'"~^^y) > is 

somehow essential given that declustering wave may appear in our model, since 
the pressure p can change from positive value to instantaneously due to a 
declustering wave. The assumption — Hme_j.o epi((p")'^) > is also needed 
in case of the appearance of new congestion regions. In regions where p""*"^ > 
and p" > 0, we have naturally that {p"^^Y — > p* and (p")'^ — > p* as e — > 0. 



Taking the divergence of (20) and using (19 1, we can indeed recover (14), which 
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leads to the incompressibility of (q""'"^)'^: • (qr"+i)"+i = 0. Here the implic- 
itness in (14) is crucial. Then ( 11 ) follows. Although we can only prove the AP 
property inside a congestion region, the numerical solutions provide evidence 
that the scheme is globally AP, including at transition between compressible 
and incompressible region. 

As for the Gauge method, it is also AP since it is a direct consequence of 
^ and ([15|. 

In the numerical simulations, we will check the AP property for concrete 
test cases. 



4 Full time and space discretization 

In this section, we present the one dimensional full time and space discretization. 
The two dimensional discretization is the easy extension of the one dimensional 
case. For the sake of completeness, we include it in the appendix. 

The Direct method: In the following, we consider the domain [b, c] — [0, 1]. 
Let the uniform spatial mesh be Aa; = jj, where M is a positive integer. Denote 



hyU: 



n+l 



the approximations of U = {p,q) at time t" — nAt and 



positions xj — jAx, for j = 0, 1, . . . , M. We fully discretize the scheme ( 14 1 and 
( 15 ) in the spirit of a local Lax-Friedrichs (or Rusanov) method [57] as follows: 



At 



1 
Ax 



Q J + 1 / 2 ( t^" , i,U"^^, Uj'J^^ ) 

/2(c/;_i,[/;,c/;_+\c/;+^) 



= 0, 



At 



Ax 



2Ax 

where the fluxes are given by 



^n+1/2 1 r n+l , ri+ll ^ /^n i n n\ 
Q, + l/2 = ^ [lj + <lj + l \ - 9C^j + l/2(Pj + l - ), 



^j + 1/2 - 2 



(9?)' , [ll+l) , / n \ . / n\ 

+ + £Po{p]+i) + epoip]) 

Pj P3+1 



(34) 



(35) 



(36) 



(37) 



They consist of the sum of a centered flux, implicit in (36 1 , explicit in (37) and 
biased terms introducing diffusion. The quantity (^"+1/2 ^^'^ local diffusion 
coefficient and is given by: 



c: 



i+1/2 
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p7 



(38) 
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It is defined as the local maximal characteristic speed related to the explicit 
pressure pq. Therefore, it remains bounded as e goes to zero. The quantity 



C 



J + l/2 



provides a numerical viscosity which is needed for scheme stability. We 
note that only the central discretization part of the flux is taken implicit, while 
the numerical viscosity part is kept explicit. In the momentum flux, only the 
part of the flux which relates to pressure is taken implicit. 

Based on this discretization, we can apply the same strategy as described in 
section 



3.2 to get an elliptic equation in p by substituting (35) into (34 1 
p1 - At 



At 



2Ax 

1 



4Aa;2 



2Ax 
At 
2Ax2 



[Cj + l/2{Pj + l - P'j ) - Cj-l/2{Pj - Pi-l)] 



(39) 



pn _ pn _ pn 



f: 



J -3/2 



= 0. 



This equation is consistent with equation ( 19 ) of the Direct method. Then we 
get a nonlinear equation for pi : 



^,Hp^)lU-nPiyr+^(p^)1^2] 



2Ax ^^J+i 
At^ 



2Ax 



P^) - C,-l/2ip'} - P^-^ 



2Aa;2 



pn pn , pn 

^j+3/2 ^J + 1/2 + 



(40) 



As mentioned before, we will use Newton iterations to solve this nonlinear equa- 
tion and get Pi^^ ■ The density p""*"^ is then obtained by inverting the nonlinear 
function pi = pi{p) with another Newton iteration. Once p"'^^ is solved, q""*"^ 
can be obtained by 



= $(c/;_i,c/;,c/;Vi)- 



with 



$(c/;Li,l/;\c/;Vi) = 9" 



At 
2Ax 

At 
Ax 



[ep,{p^;tl)-^pM-l)]' 

^J + 1/2 ^j-1/2 ■ 



(41) 
(42) 



The Gauge method Also based on (34) and (35 ), we can have the full time 



space discretization of the Gauge method. Indeed, for the Gauge method, we 



need to discretize ( 40 ) . This leads to 
1 



4Aa;2 



-,"+1 



2^;+i + ^^+i] 



1 



At 



' 2Aa; 

(q" «) g") 



[C,+i/2(p"+i - P") - C,-i/2(p," - pU)] ' (43) 



+ epo(p")+epi(p"+') 
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2 ^ 

1 



1 



2A2; ^^^+ 



(44) 
(45) 



The above equations are the direct consequences of (24)- (25 1, (28) and (32). 
However, in the numerical simulation, we will mainly test the schemes with (43 ) 
replaced by 



1 



Aa;2 



— (p"+' 



Pi) 



2Ax 



[C,+i/2ipl+i - pD - C,-i/2ipl - pU)] ' 

(46) 



which may be justified as being the direct discretization of (25) with some 
numerical viscosity added to the right hand side. The stencils are different 
in two cases. We call the Gauge method with stencil (43) Gauge 2 method 
and the one with stencil (46 1 Gauge 1 method, since they use grids <fj±2 and 
respectively in addition to ipj. There is a big difference in performance 
between the two discretizations. In fact, we will see in the next section that 
the Gauge 2 method ( with (43)) yields almost the same numerical result as the 
Direct method, while the Gauge 1 method (with ([46|) performs quite differently 
from the Direct method. This may be partially due to the fact that the Gauge 
1 method introduces more diffusion than the Gauge 2 method, which can be 
seen by inserting the Taylor expansion of 'p1±i , ¥'^±2 around x = jAx into the 
discretizations: 



1 



Ax 



2 lv'j+i 



2^^+i + 



1 



4Aa;^ 



2^"+i 



idx^ 



93"+^(jAx)Aa;^ + 0(Ax* 



Numerical diflFusion An important issue about the scheme is the numerical 
diffusion. From the equation (40), it can be seen that the diffusion for p is of 
the order of 

^{\u"\ + ^ep'^{p^))Ax + Atep[{p^)\ A^p" + AtA^{p^u^ ® u"). (47) 



And similarly, by inserting ( 34 ) into the pressure pi term in ( 35 ) we can see 
that the diffusion for q is of the order of 



\{\u^\ + ^ep'^ip^))^ AxA^q" + AteK(p")A.<7". 



(48) 



To damp out the oscillations in the mass and momentum equations, the required 
numerical diffusion is |27j 

Q(|u"| + v/^^V^)) AxA.p", Q(|ii"| + V^^))^ AxA.g", (49) 
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respectively. To ensure that this numerical diffusion is achieved, we need the 
condition 



\{\u-\ + ^7^))/^x + ep',{p-)M\ > + A.T, (50) 



which leads to 




2(V^p;,(p") + v/ep'(p") 



1 



) 



(51) 



This condition is automatically satisfied in the congested region (p — >■ p*) for 
small e, since ep'{p'^) — >■ oo as e — > 0. However, it contradicts the CFL condition 
in the non-congested region for small e, since ep'[p"-) and epo(p") — as e — 0. 
From this analysis, there should be no oscillations in the congestion region 
while the numerical diffusion may not be sufficient in the non-congested region. 
However, numerical simulation seems to indicate that the numerical viscosity 
in this scheme is sufficient to damp out the oscillations in the non-congested 
region. 

5 Numerical results 

5.1 One dimensional test cases 

In this section, we use several numerical examples to test the performance of 
the schemes. Corresponding to different situations, four examples are tested. 
All these examples are the compositions of Riemann problems. Since the exact 
solutions to Riemann problem can be determined as in section[2] we can compare 
the exact and numerical solutions. Different measurements of the relative errors 
will be applied to test the performance of our schemes. And the numerical 
Courant number is computed. 

In the following, we choose 7 = 2 and the maximal density p* — 1. The test 
problems are: 




X £ (0.75,1], 
X e [0,0.5), 

X e (0.5,1], 



X e [0,0.5), 
X e (0.5,1], 



X e [0,0.5), 
X e (0.5,1], 

X e [0,0.25), 



X e (0.25,0.75), 



(55) 



(53) 



(54) 



(52) 
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0.5 



-0.5 





— — - Direct 




>»« — 


Exact 
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X 





0.2 0.4 0.6 0.. 



Figure 3: The direct and Gauge 1 methods for problem (PI) at t = 0.05 with 
£ = 10~*, Ax = 5 X 10~^, At = 5 X 10^*. The soUd hues are the exact solutions. 
The dashed curves are the numerical solutions. The left graphs are for p, the 
right ones for q, both as functions of x. 



The first example (PI) illustrates how the AP schemes capture shocks near 
congestion. The second example (P2) shows how the AP schemes work near 
vacuum. The third example (P3) simulates the interaction of two shocks near 
congestion. The last example (P4) shows some problems in the Gauge 1 method 
and will be used to justify the splitting of p. 

Example 1. The solution to the Riemann problem (PI) consists of two 
shocks propagating in the opposite directions. The density of the intermediate 
state is close to the maximal density. In the following, we will test the two 
methods described in section [4] with different parameters e and different mesh 
sizes Ax, At. 

1. First, we choose e = 10-*,Aa; = 5 x 10'^, At = 5 x 10"''. We wiU 
compare the performances of the two methods proposed in section [4] It 
can be seen from Figure [3] that there is large oscillations of the momentum 
in the congested region. But the propagation of the shock is captured well. 
In comparison, the Gauge 1 method as illustrated in Figure [3] eliminates 
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0.2 
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Figure 4: Fix Aa; = 5 x 10 At = 5 x 10 *. The numerical results of the 
Direct method for problem (PI) aX t = 0.05 with e = 10"^ and e = 10~*. The 
left graphs are for p, the right ones for g, both as functions of x. 



all the oscillation. 

2. We look how the choices of the parameters e Ax and At affect the numer- 
ical result. We may fix Ax, At but change the value of e so as to test the 
cases e < At and e > At. From the numerical results, it can be seen that 
the oscillations in the momentum always appear for different choices of e 
but are smaller as e — ^ for this choice of parameter. This verifies the AP 
property. As for the Gauge 1 method, it has the same performance for all 
value of £. Thus it shares the same property. 

3. The above observation can be quantitatively investigated by measuring 
the difference between the numerical solution W and the theoretical one 
w. We use two measurements: one is the relative error of the numerical 
solution W compared with w in the sense of L"'^ norm and the other is the 
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Parameters 



Direct 



Gauge 1 



£ 


Ax 


At 


e{W) 




ratio 


e{W) 




ratio 




1/200 


1/250 


1.1012 X 10" 


-2 




1.1209 X 10" 


-2 






1/200 


1/500 


8.0103 X 10" 


-3 


1.3747 


7.8670 X 10" 


-3 


1.4248 




1/200 


1/1000 


4.2631 X 10" 


-3 


1.8790 


4.9486 X 10" 


-3 


1.5897 




1/200 


1/2000 


5.1528 X 10" 


-3 


0.8273 


5.4107 X 10" 


-3 


0.9146 


10"'' 


1/200 


1/10000 


5.6843 X 10" 


-3 


0.9065 


1.5761 X 10" 


-2 


0.3433 


1/200 


1/1000 


4.2631 X 10" 


-3 




4.9486 X 10" 


-3 






1/400 


1/2000 


3.5612 X 10" 


-3 


1.1971 


3.1706 X 10" 


-3 


1.5608 




1/800 


1/4000 


1.3085 X 10" 


-3 


2.7216 


1.4253 X 10" 


-3 


2.2245 




1/1600 


1/8000 


5.7676 X 10" 


-4 


2.2687 


6.1850 X 10" 


-4 


2.3044 




1/1600 


1/16000 


7.2302 X 10" 


-4 


0.79771 


7.0713 X 10" 


-4 


0.87466 




1/200 


1/250 


1.3188 X 10" 


-2 




1.3721 X 10" 


-2 






1/200 


1/500 


7.7490 X 10" 


-3 


1.7019 


8.5433 X 10" 


-3 


1.6061 




1/200 


1/1000 


7.6793 X 10" 


-3 


1.0091 


7.8588 X 10" 


-3 


1.0871 




1/200 


1/2000 


8.2751 X 10" 


-3 


0.9280 


7.5447 X 10" 


-3 


1.0416 




1/200 


1/10000 


1.0758 X 10" 


-2 


0.7692 


1.5761 X 10" 


-2 


0.4787 


1/200 


1/1000 


7.6793 X 10" 


-3 




7.8588 X 10" 


-3 






1/400 


1/2000 


3.7602 X 10" 


-3 


2.0423 


3.8637 X 10" 


-3 


2.0340 




1/800 


1/4000 


1.7934 X 10" 


-3 


2.0967 


1.8646 X 10" 


-3 


2.0721 




1/1600 


1/8000 


8.0483 X 10" 


-4 


2.2283 


8.7074 X 10" 


-4 


2.1414 




1/1600 


1/16000 


7.9308 X 10" 


-4 


1.0148 


8.4647 X 10" 


-4 


1.0287 



Table 1: Comparison of the relative error between the Direct and Gauge 
1 methods at t = 0.025. The 'ratio' column ratio provides comparisons of the 
relative norm error between the previous and current rows, where either Ace 
is fixed or Ax/ At is fixed. 



difference of their total variation: 

e{W) = where |1^|U. = ^El^^^l' (^6) 

j 

|Tot.Var.{W^}-Tot.Var.{w}| , rx. r i V^, 

g{W) = J Tot VarW ' Tot. Var.{z«} = E|u;,+i - 

i 

(57) 

In Table [l] and [2j the relative error in terms of distance and total 
variation between the numerical and exact solutions at < = 0.025 are 
computed. The Gauge 2 method yields quite similar result to the Direct 
method. So it is not listed in the table. It can be seen that the Direct 
method is usually better than the Gauge 1 method in the norm. To 
reflect the observation we made from looking at Figure [3j we use the 
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Parameters 



Direct 



Gauge 1 



s 


Ax 


At 


9{W) 




ratio 






ratio 




1/200 


1/250 


7.4909 X 10" 


-3 


- 


7.2225 X 10" 


-4 


- 




1/200 


1/500 


7.5636 X 10" 


-2 


0.0990 


6.7956 X 10" 


-3 


0.1063 




1/200 


1/1000 


4.1433 X 10" 


-1 


0.1826 


2.4506 X 10" 


-2 


0.2773 




1/200 


1/2000 


8.5937 X 10- 


-1 


0.4821 


1.4766 X 10" 


-2 


1.6596 


10-^ 


1/200 


1/10000 


1.0603 




0.8105 


1.8995 X 10- 


-2 


0.7774 


1/200 


1/1000 


4.1433 X 10" 


-1 




2.4506 X 10" 


-2 






1/400 


1/2000 


4.0126 X 10" 


-1 


1.0326 


3.8585 X 10- 


-2 


0.6351 




1/800 


1/4000 


2.5996 X 10" 


-1 


1.5435 


2.4447 X 10- 


-2 


1.5783 




1/1600 


1/8000 


5.4191 X 10" 


-1 


0.4797 


2.3520 X 10- 


-2 


1.0394 




1/1600 


1/16000 


1.0294 




0.5264 


1.3054 X 10- 


-2 


1.8017 




1/200 


1/250 


5.1377 X 10- 


-5 




4.1271 X 10- 


-4 






1/200 


1/500 


5.2584 X 10" 


-3 


0.0098 


2.1808 X 10- 


-4 


1.8925 




1/200 


1/1000 


1.1354 X 10^ 


-1 


0.0463 


2.8711 X 10- 


-3 


0.0760 




1/200 


1/2000 


4.9598 X 10" 


-1 


0.2289 


9.9861 X 10- 


-3 


0.2875 


10-2 


1/200 


1/10000 


1.2766 




0.3885 


4.4336 X 10- 


-3 


2.2524 


1/200 


1/1000 


1.1354 X 10" 


-1 




2.8711 X 10- 


-3 






1/400 


1/2000 


1.1856 X 10" 


-1 


0.9577 


2.9617 X 10- 


-3 


0.9694 




1/800 


1/4000 


1.2408 X 10" 


-1 


0.9555 


2.9916 X 10- 


-3 


0.9900 




1/1600 


1/8000 


1.2200 X 10" 


-1 


1.0170 


2.7207 X 10- 


-3 


1.0996 




1/1600 


1/16000 


3.9234 X 10" 


-1 


0.3110 


1.3249 X 10- 


-2 


0.2054 



Table 2: Comparison of the total variation relative error between the Direct 

and Gauge 1 methods at i = 0.025. The 'ratio' column provides comparisons of 
the total variation relative error g{W) between the previous and current rows, 
where either Aa; is fixed or Ax /At is fixed. 
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e 




Ax 


stable At 


max A 


(maxA)ll 


10- 


-4 


1/100 


1/100 


1.2186 


1.2186 


10- 


-4 


1/200 


1/210 


1.2669 


1.2066 


10- 


-4 


1/400 


1/430 


1.1962 


1.1128 


10- 


-4 


1/800 


1/870 


1.1962 


1.1000 


10- 


-4 


1/1600 


1/1730 


1.1938 


1.1041 


10- 
10- 
10- 
10- 
10- 


-2 
-2 
-2 
-2 
-2 


1/100 
1/200 
1/400 
1/800 
1/1600 


1/100 
1/280 
1/580 
1/1170 
1/2340 


1.7903 
1.6571 
1.6486 
1.6486 
1.6645 


1.4902 
1.1836 
1.1370 
1.1272 
1.1381 



Table 3: The numerical Courant number for the Direct method, max A is the 
maximal eigenvalue of the explicit part of the scheme for all time steps before 
the waves reach the boundary. 



total variation norm, which captures the oscillations better. For the total 
variation measurement, it can be seen that the Gauge 1 method is always 
better in controlling the oscillations. Another observation can be made 
from the tables is how the accuracy is changed with different parameters. 
For both two measurements, we test the accuracy with Ax fixed or Ax/ At 
fixed. In the test where Ax is fixed, it can also be seen that the relative 
error in norm can not be reduced much by refining the time mesh from 
1/500 to 1/10000 with fixed Ax = 1/200. This feature is the same as 
the standard hyperbolic solvers: better accuracy can not be obtained by 
using a smaller At once the scheme is stable. In the test where Ax/ At 
is fixed, we check how the relative error is decreasing with respect to 
Ax. It is interesting to see that the error is not always decreasing. And 
since the convergence order for explicit local Lax-Friedrichs scheme with 
discontinuities is ^ [17] , we may expect that the relative error in norm 
is reduced by -\/2 when the space mesh is refined by 2 with Ax/ At fixed. 
However, this is not the case. 

4. We will also check the numerical Courant number in tables [3] and ID It can 
be seen that the CFL condition of our scheme is greatly improved from 
the one for standard hyperbolic solver At = 0(£Ax). 

5. We can also quantify the observation as in Figure |4| In Table[5] the relat ive 
errors of solutions in the norm or total variation are summarized. The 
data for e = 10"^ lO"'' are the same as those in Table [l] and [2] It can 
be seen that there is no big increase in the relative error for different e 
with fixed Ax, At. As discussed in section]^ the theoretical solutions to 



system (62) and (63) with positive e are converging to the solutions to 
systems with e = when e — > 0. Thus, the numerical solutions tend to 
the theoretical solutions systems with e = for fixed Ax, At as well. This 
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e 




Ax 


stable Ar 


max A 


max A^c— 


10- 


A 
-4 


1/100 


1/100 


1.2151 


1.2151 


10- 


-4 


1/200 


1/210 


1.2279 


1.1695 


10- 


-4 


1/400 


1/500 


1.3482 


1.0786 


10- 


-4 


1/800 


1/1340 


1.4398 


0.8596 


10- 


-4 


1/1600 


1/3130 


1.3555 


0.69293 


10- 


-2 


1/100 


1/100 


1.6480 


1.6480 


10- 


-2 


1/200 


1/260 


1.6475 


1.2673 


10- 


-2 


1/400 


1/540 


1.6475 


1.2204 


10- 


-2 


1/800 


1/1190 


1.6530 


1.1112 


10- 


2 


1/1600 


1/2520 


1.6490 


1.0470 



Table 4: The numerical Courant number for the Gauge 1 method, max A is the 
maximal eigenvalue of the explicit part of the scheme for all time steps before 
the waves reach the boundary. 



versifies the AP property. 

6. With a slightly variant version of this test case, we can also see the role 
played by the splitting of pressure p. By making the momentum of the left 
and right states 10 times smaller, we have another test case (PI') sharing 
the similar behaviour of (PI). 

In Figure [5j the numerical solutions obtained by the Direct method with 
and without po are compared. This confirms the observation made in 
[l8] for the low Mach number limit. It can be seen that a lot of extra 
oscillations appear in the numerical solutions when there is no splitting of 
pressure p (po = 0,pi = p). That is the reason why we should add a pq 
term in the explicit part of the numerical scheme. 

Example 2. The Riemann problem (P2) is obtained by exchanging the left 
and right states of Riemann problem (PI). So the solution to the problem (P2) 
consists of two rarefaction waves and a vacuum state appears as the intermediate 
state. As shown in proposition |4j these two rarefaction waves tend to be contact 
waves. 

As in example 1, We choose e = 10~'',Aa; = 5 x 10~^,At = 5 x 10~^. 
t=0.0005. It can be seen from Figure [6] that the Direct method captures the 
vacuum and rarefaction waves well. In comparison, the Gauge 1 method as 
illustrated in Figure [6] shows larger diffusion. 

Example 3. The solution to the problem (P3) consists of two Riemann 
problems: both of them are like the Riemann problem in (PI). So there are two 
congested regions and eventually they will collide. We are interested in observing 
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Figure 5: The numerical results of the Direct method with and without po for 
problem (PI') at t = 0.2 with e = 10"'', Ax = 5x10"^ At = 5x10"''. The solid 
lines are the exact solutions. The dashed curves are the numerical solutions. 
The left graphs are for p, the right ones for q, both as functions of x. 
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Figure 6: The Direct and Gauge 1 methods for problem (P2) aX t = 0.05 with 
e — 10"'*, Ax — 5x 10~^, At = 5 x 10""*. The soUd Unes are the exact solutions. 
The dashed curves are the numerical solutions. The left graphs are for p, the 
right ones for q, both as functions of x. 
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Parameters 



Direct 



Gauge 1 



£ 


Ax 


At 


e{W) 




g{w) 




e{W) 




9{W) 






1/200 


1/1000 


7.6793 X 10- 


-3 


1.1354 X 10" 


-1 


7.8588 X 10- 


-3 


2.8711 X 10" 


-3 


10-2 


1/200 


1/2000 


8.2751 X 10- 


-3 


4.9598 X 10- 


-1 


7.5447 X 10" 


-3 


9.9861 X 10" 


-3 




1/800 


1/4000 


1.7934 X 10- 


-3 


1.2408 X 10- 


-1 


1.8646 X 10- 


-3 


2.9916 X 10- 


-3 




1/1600 


1/16000 


7.9308 X 10- 


-4 


3.9234 X 10- 


-1 


8.4647 X 10- 


-4 


1.3249 X 10- 


-2 




1/200 


1/1000 


4.2631 X 10- 


-3 


4.1433 X 10- 


'1 


4.9486 X 10- 


-3 


2.4506 X 10" 


-2 


lo-'' 


1/200 


1/2000 


5.1528 X 10- 


-3 


8.5937 X 10- 


-1 


5.4107 X 10- 


"3 


1.4766 X 10- 


-2 




1/800 


1/4000 


1.3085 X 10- 


-3 


2.5996 X 10- 


-1 


1.4253 X 10- 


-3 


2.4447 X 10- 


-2 




1/1600 


1/16000 


7.2302 X 10- 


-4 


1.0294 




7.0713 X 10- 


-4 


1.3054 X 10" 


-2 




1/200 


1/1000 


7.0600 X 10- 


-3 


6.4521 X 10- 


-1 


5.8159 X 10- 


-3 


1.4861 X 10" 


-2 


10-8 


1/200 


1/2000 


5.8872 X 10- 


-3 


6.4091 X 10- 


-3 


5.8872 X 10- 


-3 


6.4091 X 10- 


-3 




1/800 


1/4000 


1.7611 X 10- 


-3 


6.7218 X 10- 


-1 


1.4624 X 10- 


-3 


1.5544 X 10- 


-2 




1/1600 


1/16000 


9.1296 X 10- 


-4 


8.2410 X 10- 


-1 


8.0587 X 10- 


-4 


1.6092 X 10- 


-2 



Table 5: The and total variation relative error of the Direct and Gauge 1 
methods at i = 0.025 for different e. 



how the numerical methods behave at coUision. We fix Ax — bx 10^"^, At = 
5 X 10-^ and choose e = 10~^ and 10-^. Since only shocks are involved, we will 
use the Gauge 1 method only. 

From Figures [7] and [8j as e becomes smaller, it takes shorter time to form a 
new congestion region from the two colliding congestion regions; from 48 time 
steps to no more than one time step. It can be seen that as £ — > 0, the collision 
of these two congested shocks aggregate instantaneously. 

Example 4. The solution to the Riemann problem (P4) consists of two 
shocks with intermediate state away from the congestion density. So in the 
second Riemann problem there are two rarefaction waves and a vacuum state 
appears as the intermediate state. As shown in section [2] these two rarefaction 
waves tend to be contact waves. 

As above, we choose e = 10^^, Ax = 5 x 10^'^,At = 5 x 10^^. It can be 
seen from Figure |9] that the Direct method performs well when the density is 
far away from congestion. However, the Gauge 1 method does not work well 
possibly due to the extra diffusion. 

Remark 2 All the features described above are preserved when we apply the 



fully implicit method (21](22) (implicit in both p and to the Direct and 



Gauge 1 method. There is a little improvement in the accuracy but no major 
one. 

Finally, we can compare the Gauge 1 and Gauge 2 methods. With the stencil 



of (46) in the same setting as that of example 2, the Gauge 2 method yields 



almost the same numerical result as the Direct method (see Fig. 10 ) 
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Figure 7: Gauge 1 method for problem (P3) with e = 10 Ax = 5 x 10 ^,At = 
5 X 10^^. The left graphs are for p, the right ones for q, both as functions of x. 
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Figure 8: 

5 X 10-4 



Gauge 1 method for problem (P3) with e = 10 ^, Aa; 5 x 10 ^,At = 
The left graphs are for p, the right ones for q, both as functions of x. 



26 










1 


1 




/ 


1 




. / 


1 






1 






\ 






\ 








\ 






\ 






\ 









0.2 0.4 0.6 0.8 

X 



0.35 
0.3 

0.25 
0.2 

0.15 
0.1 





— — — Gaugel 

Exact 


s 




\ 




\ 




( 




\ 


I 

\ 

\ 

\ 







0.2 0.4 0.6 0.8 1 

X 



Figure 9: The Direct and Gauge 1 methods for problem (P4) at t = 0.2 with 
e = 10"'', Ax = 5 X 10-3, At = 5 x lO""*. The sohd hues are the exact solutions. 
The dashed curves are the numerical solutions. The left graphs are for p, the 
right ones for q, both as functions of x. 
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Figure 10: The Direct and Gauge 2 methods with discretization (43 ) for problem 
(P2) at t — 0.05. The solid lines are the exact solutions. The dashed curves are 
the numerical solutions. The left graphs are for p, the right ones for q, both as 
functions of x. 
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5.2 A two dimensional test case 



In this section, we will test the Direct and Gauge 1 method in the 2D case. 
Since there is no theoretical solution in 2D case, only some phenomena will be 
presented. 

The test example is chosen to illustrate the collision of two congested regions. 
It is basically a two dimensional extension of the test case (P3) with some lateral 
shift. The initial data of the density and velocity is: 



p = 0.8 X IauB + 0.6 X l[o,l]x[0,l]\(yluB)i 







A = 


'1 5 ■ 

6' 12 


X 


'1 7 " 

3' 12 



B 



-1 


7 5 
12' 6 



1; 





■ 5 2" 


X 


12' 3 



(59) 
(60) 

(61) 



The vector field q is plotted in the Figure [TT] 



Remark 3 In this test example, the background density is p = 0.6. It is even 
more interesting to see what happens when the background density is close to 
zero. However, our schemes are not performing well in this situation. Indeed, 
vacuum is a big problem which needs some special treatment. And when the 
background density is close to zero, the system is almost pressureless, which is 
also a difficult problem. These problems may be considered in the future work. 

It can be expected that there will be two congested regions forming and 
moving towards each other with two shocks in the front and two rarefaction 
waves left behind. These two shocks compress the fluid and cause congestion. 
This is reflected in Figure [T2j 

These two shocks with opposite directions meet at a later time. The interac- 
tion of these shocks forms a bigger congestion region with higher density, which 



can be illustrated in Figure 13 The interaction of the two shocks induces the 
formation of two new shocks moving in the orthogonal direction compared to 
the initial motion. 

The similar result can be obtained by the Gauge 1 method as in Figure [T4| 
and[l5l 

The difference between the two methods in the two dimensional case is not 
striking. In order to illustrate them, we look at a cut at y = 0.5. The cuts are 
displayed on Figure [16] and [TT] 

Similar observations as in the one dimensional case can be made from Fig- 
ure [16] and [17] The Gauge 1 method provides a little less oscillations in the 
congestion region but also brings more diffusion. 
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Figure 11: The initial data of the density and momentum q, both as functions 
of X and y. 



30 



0.2 



0.4 



0.6 



0.8 



(a) density p 

1p- — 

0.8-: 



0.6 
0.4 



0.2 



0.2 0.4 0.6 0.8 1 

X 

(b) vector field of momentum q 

Figure 12: The Direct method with e = 10""*, Ax = 5 x 10"^, At = 5 x 10""* 
at t = 0.05. The left graph is for the right one for vector field g, both as 
functions of x and y. 
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(b) vector field of momentum q 

Figure 13: The Direct method with e = 10"*, Ax = 5 x lO"^^ A< = 5 x 10""* at 
t — 0.2. The left graph is for p, the right one for vector field q, both as functions 
of X and y. 
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Figure 14: The Gauge 1 method with e = 10""*, Aa; = 5 x 10"^ At = 5 x 10""* 
at t = 0.05. The left graph is for p, the right one for vector field q, both as 
functions of x and y. 
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(b) vector field of momentum q 



Figure 15: The Gauge 1 method with e = 10""*, Aa; = 5 x 10"^ At = 5 x IQ-^ 
at t = 0.2. The left graph is for p, the right one for vector field q, both as 
functions of x and y. 
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Figure 17: Momentum of the Direct and Gauge 1 methods with e = 10~^, Ax 
5x 10"'^, At = 5x 10~^ att = 0.05. Cut of the momentum along the hue y = 0. 
as a function of x. Left hand: Direct method; right hand: Gauge 1 method. 
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6 Conclusion 



In this paper, we have studied the Euler system with a maximal density con- 
straint. A small parameter e was introduced to measure the stiffness of the 
constraint. As e — >■ 0, the model gives rise to a two-phase model with congested 
regions (with maximal density) and uncongested regions (with density below 
the maximal density). 

One-dimensional solutions of this asymptotic problem have been investigated 
to provide the information of the interface conditions. However, it can not char- 
acterize the whole dynamics and it is hard to extend to higher dimensional cases. 
Therefore, we have devised asymptotic preserving numerical schemes, which are 
valid for all range of e and thus are capable of capturing the asymptotic dy- 
namics. Two numerical schemes have been considered and compared on both 
one-dimensional and two dimensional test-cases. They both capture the con- 
gested regions well. However, the first method shows some oscillations near 
the interface between congested and uncongested regions, while the second has 
much less oscillations but is more diffusive. A careful error analysis in differ- 
ent norms with respect to different parameters (time and space steps, e) are 
conducted for a one dimensional test case. 

Following [35, , the development of second order schemes will be performed in 
future works. More robust schemes for capturing the divergence free constraint 
in the congested regions could also be investigated. Coupling our methodology 
with schemes dedicated to the pressureless gas dynamics could also improve the 
results in the low density regions. Finally, simulations of the non-conservative 
model with a supplementary geometric constraint on the speed of the flow |17) 
will be an interesting problem. 
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CT-2005-021122), by the Agence Nationale de la Recherche (ANR) under con- 
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A Appendix : Solutions of the one-dimensional 
problem: the Riemann problem and the clus- 
ter collisions 

The one-dimensional version of system Q-© can be written as follows: 



where q{x, t) is here a scalar function and x is the position in M. In this section, 
we first investigate the Riemann problem and the limits of its solutions as e — > 



dtp + d^q 



= 0, 



(62) 




(63) 
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and then wc briefly recall the one-dimensional solutions of the limit system ([s])- 
([9|-([l0l) provided in which consists of the collision of two finite clusters 
(domain where p = p*). 

A.l The one-dimensional Riemann problem 

The Riemann problem is an initial value problem, where the initial condition is 
a piece-wise constant function with a discontinuity between two constant states: 



{pQ,qo){x) = 



{pi,qi), for a; < 0, 
{pr,qr), for a; > 0. 



The solutions of this problem for the Euler system (62)-(63l are well known. 
So, the strategy is to take the limits of these solutions as e goes to zero. This 
provides the solutions of the Riemann problem for the singular asymptotic limit 
([8|-([9|-([l0]). Similar studies was carried out in [S] for the isentropic Euler equa- 
tion, in [9 for a traffic jam model and in |17) for a herding problem. 

A. 1.1 Shock and rarefaction waves 

The material of this section is classical, and is given here only for the reader's 
convenience. 

The two characteristic speeds A§_ and the two characteristic fields r!|_ of the 



one-dimcnsional system (62|-(63) are 



It can be easily checked that both characteristic fields are genuinely non linear 
for positive densities. Therefore, the solutions of the Riemann problem are made 
of constant states separated by rarefaction or shock waves |27| . 

The rarefaction waves are continuous self similar solutions: (p{x/t), q{x/t)). 
Given a state (p, q) , the states which can be connected to {p, q) by a rarefaction 
wave are those located on the integral curves i^. of the right eigenvectors of the 
Jacobian matrix of the flux function issued from (p, q) . They are given by: 

= = u± Vep'(p), p(0)=/5, i|(0) = <7, 

which is equivalent to (p) — u± ^ep'(p), i^. {p) = q and then to equation: 
z| (p) ^pu± pV^{P{p) -P{p)), (64) 



where P is an antiderivative of \Jp'{u)/u and u = q/p. The graph of (resp. 
i\) is called the 1-integral curve (resp. the 2- integral curve). The following 
proposition provides several features of the integral curves. 

Proposition 1 (Integral and rarefaction curves) 
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1. The 1-integral curve i'L (resp. the 2-integral curve i^) is concave (resp. 
convex) as functions of p and «1(0) = i+(0) = 0. 

2. The limit of the integral curves as e goes to zero is the union of the straight 
lines {q = pii} and {p = p*}. 

3. Suppose that the state (f is such that (f — > p* and epijf) — >■ p. For all 
p < , we have: 



u 



<V-er ^du= O {e^) 
Jo u P'^P' 



4- If{p, q) is a left state, the right states which can be connected to it by a rar- 
efaction wave are those located on the 1-rarefaction curve {{p, iL{p)), p < p} 
or the 2-rarefaction curve {{p,i^{p)), p > p}. 

The last point of this proposition stems from the compatibility conditions of the 
characteristic speeds. The proof of this proposition is classical and omitted. 

A shock wave is a discontinuity between two constant states, {p, q) and 
(p, g), travelling at constant speed a. The states {p. q), which can be connected 
to (p, q) by a shock wave, are determined by the Rankine-Hugoniot conditions: 



^ + ep(p) 



a[q], (65) 



where [/] := f — f for all quantities /. Easy computations show that the 
admissible states are of the form (p, /i§_(p)), where is : 

h% {p)=u± ^^^/{p-p){ep{p)-Ep{p)). (66) 

and the shock speeds are: 




{p- p) V p V (p- p) 



The graph of h'L (resp. /i^) is called the 1-Hugoniot curve (resp. the 2-Hugoniot 
curve). The following proposition provides several properties of the Hugoniot 
curves: 

Proposition 2 (Hugoniot and shock curves) 

1. The 1- Hugoniot function ht_ (resp. the 2-Hugoniot function h^) is concave 
(resp. convex) and /il(0) = ^.+ (0) = 0. 

2. The limits of the graphs of lf__ and h^_^_ when e goes to zero are the union 
of the straight lines {q = pii} and {p = p*}. This is true also when p = p^ 
depends on e and tends to p* . 
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Figure 18: Curves issued from {p,q) = (0.5,0.5) and s £ {lO, 1,10~^}. In 
dotted line : the rarefaction curves. In solid line : the shock curves. In blue : 
the 1-curves. In black : the 2-curves. As e — ^ 0, the curves are tending to the 
straight lines {q = up} and {p — p*}. Parameters : k = 2, p* — 1. 



3. If (p, q) is a left state, the right states which can he connected to it by an en- 
tropic shock wave are those located on the 1-shock curve {(p, h'L{p)), p > p} 
or the 2-shock curve {(p, /i'^(p)), p < p}. 

The proof of this proposition is classical and omitted. 



A. 1.2 The Riemann problem for system (62|-(63) 



Geometric considerations on the intersection of the integral and Hugoniot curves 
enable us to solve the Riemann problem j27| . These arguments are really sim- 
plified in the limit £ — ?> 0, due to the much simpler behaviour of the integral and 
Hugoniot curves: they both converge to union of the straight lines {q — pu\ and 



{p — p*}. This behaviour is illustrated in Fig. 18 The following proposition 
provides the nature of the solutions of the Riemann problem for e small enough. 
It depends on the sign of the relative velocity ui — Ur, where ui = qe/pe and 

Ur = qr/Pr- 

Proposition 3 Let {pi,qi), (Pr^qr): ^^ft o,i^d right states. Then the solutions 



of the Riemann problem related to 1(6^- {63) have the following forms 



1. If Ui < Ur, then for e small enough, the intermediate state is the intersec- 
tion point of the 1-integral curve issued from (pi,qi,pi) and the 2-integral 
curve issued from {pr,qnPr)- Besides, the intermediate density p is lower 
than pr and pg. This is valid even if pi or/and pr tends to p* . 
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2. If ui > Ur, then for s small enough, there are two subcases: 

• if h'L{Pr) > Qr md < Qi, the intermediate state is the in- 
tersection point of the 1-Hugoniot curve issued from [p^, qg) and the 
2-Hugoniot curve issued from {pr,qr)- 

• */ h'L{Pr) < Qr (resp. h^-^{p\) > qi), the intermediate state is the 
intersection point of the 1-Hugoniot curve (resp. 1-integral curve) 
issued from {p^,q^) and the 2-integral curve (resp. 2-Hugoniot curve) 
issued from {pr,qr). 

3. If Ur ~ U£ and pi < pr (resp. pi > pr), then the intermediate state is the 
intersection point of the 1-Hugoniot curve (resp. 1-integral curve) issued 
from {p£,q£,pi) and the 2-integral curve (resp. 2-Hugoniot curve) issued 
from [pr, qr,pr). Besides, the intermediate density p is the interval [pi, pr] 
(resp. [pr,pe]). This is valid even if pe or/and pr tends to p* . 



This proposition results from propositions [T] and [2j The proof is omitted here (it 
is an easy adaptation of the proof of Theorem 1 in 17J). Note that the nature of 
the curves (integral or Hugoniot curve) implies the nature of the waves involved 
in the Riemann problem. 



A. 1.3 The limits of solutions of the Riemann problem. 

We introduce the following initial conditions: 

IP0,90JW-| (^s^^e)^ forx>0. 

with {p\,ql,ep{p\)) {pi,qe,P£) and {Pr,qr,£p{Pr)) ^ {Pr,qr,Pr) as e goes to 
zero. 

The following proposition provides the solution when pi,pr < p* and so 
Pe = Pr — 0. The nature of the solution only depends on the sign of the relative 
velocity ue — Ur, where ue = qij pt and — qrj Pr- 

Proposition 4 (case pi < p* , pr < p* ) 

1. If Ue — Ur < 0, then the solution consists of two contact waves connecting 
the two states to the vacuum. This is summarised in the following diagram: 

, „x contact ,^ „s vacuum „s contact , „s 

{pi,qi,\)) — > (0,g^0) — > (0,gr,0) — > [pr,qr,Oj. 



2. If Ui — Ur > 0, then the solution consists of two shock waves connecting 
the left state {pi, qi, 0) to an intermediate state {p* ,q,p) and then (p*, q,p) 
to the right state (pr,9r,0); 

/ ^\ shock / ~ _\ shock , 
[Pi,qe,0) — > (p ,q,p) — > (pr,qr,0), 
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where the intermediate momentum q and the intermediate pressure p are: 

q = Uip* - \ —\/{p* - Pt)P = UrP* + \ —\/{p* - Pr)P, 
\ Pi \ Pr 



- I x2 / P - Pi , P - Pt 

VV P(P V PrP* 

and the shock speeds a- and are given by: 



-2 



^ue~ . tVp, a+=Ur . , ^ . 

V Pl[P - Pe) V Pr[P* - Pr) 

3. If ui = Ur, then the solution consists of only one contact wave connecting 
{pi,qi,0) to {pr,qr,0): 

, „^ contact , 

(P£,9£,0) > (Pr,qr,0). 

It is the same results as those obtained in [T^] , where they directly proved it by 
defining a notion of entropy solutions for the asymptotic problem. 

In all the following proofs, {p£, qi), {p, q), (pr, qr) will respectively denote the 
left, the intermediate and the right states involved in each different Riemann 
problems. Al (resp. A^) will implicitly refer to the first (resp. second) char- 
acteristic speed of the left state (resp. right state). The characteristic speeds 
related to the intermediate state are denoted: Al, A*^. In the following, and 
/il (resp. and /if^) will refer to the 1-curves (resp. 2-curves) issued from the 
left state (resp. right state). The notation [f]i = f — fg (resp. [f]r = f — fr) 
will denote the difference between the intermediate and the left (resp. the right) 
values of any quantity /. 

Proof 1. From propositionjSj the solution for small e consists of two rarefaction 
waves. The intermediate density solves equation il(p) = *+(p), that is: 

pUr ± p\/e[P[p)]r = pui ± p\/e[P{p)\t 

Since p is lower than pi and pr (and then [P[p)\r and [P{p)\f are bounded) and 
Ur — ui is not zero, it is easy to deduce that the limit solution of this equation 
is p = 0, which defines a vacuum state. Besides, Al and Af_ tends to (resp. 
A^ and A^ tends to Ur). Therefore, the limit waves are contact waves. 

2. From proposition |3j the solution for small e consists of two shock waves 
(since \\mh^_^{pi) = piUr and limL h^_{pt) = Prqi)- We have /il ^(p) — rip)i 
that is: 



pue - 
which yields 



— VIpUIMpJU = PUr + \ —\/[p\r[ep{p)]r, 
Pi V Pr 
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Since ug — Ur is different from zero, the limit intermediate pressure, limep(p'^) = 
lim[£p(/9)]£ — lim[£p(p)]r, is not zero. Thus p — >■ p* as e — >■ 0. Finally, the limit 
values of the two shock speeds can be easily inferred from equation (|67|. 

3. Supposing that pe < pr, then the intermediate density satisfies pe < p"^ < 
Pr and consequently Al , Aj_ , — Ui = Ur as e goes to zero, which yields a 
unique contact wave. 

We now consider the case where the left state is a congested state: pe = 
p* , Pr < p* and pi < +00. By a symmetry argument, the case pi < p* , pr = p* 
can be easily deduced. The limits of rarefactions waves when one state tends to 
congestion lead to the so-called dcclustering waves: 

Definition 1 ( Declustering waves ) A declustering wave consists of a shock 
wave between two congested state, with infinite speed, and with a zero pressure 
for positive time. 

The following proposition states the solutions of the Riemann problem. 

Proposition 5 (case pi ^ p* , p^ < p* ) 

1. Ifuf ^Ur < 0, then the solution consists of one declustering wave connect- 
ing the left state {p* ,qi,pi) to a congested and pressureless state {p* ,qi,0) 
and then a contact wave connecting {p* , g^, 0) to vacuum and another con- 
tact wave connecting vacuum to the right state {pr,qr,0): 

/ « _ s declust. , f contact „, vacuum contact , 

(p ,<li,Vi) — > (P ,Qi,0) — > (0,qe,0) — > (0,gr,0) — > (pr,gr,0). 

2. If ue ^ Ur > 0, then the solution consists of two shock waves connecting 
the left state {p*,qi,pe) to an intermediate congested state {p* ,q,p) and 
then connecting this intermediate state to the right state speeds {pr,qr,0)-' 

/it — \ shock / ~ _\ shock / .-.n 

\p ,qi,Pe) — > [p ,q,p) — > [Pr,qr,0), 

where the intermediate momentum q and the intermediate pressure p are 
given by: 

P*Pr / ^2 

q = pue, P=— (Ui-Ur), 

P - Pr 

and the two speed and are: 



3. If Ui — Ur, then the solution consists of one declustering wave connecting 
the left state (p* , qe,pi) to the intermediate state {p* , qe, 0) and one contact 
wave connecting the intermediate state to the right state (pr,qr,0): 

, * _ X declust. , * _x contact , 

\P ,qe,Pe) — > [P ,qe,0) — > {Pr,qr,0). 



43 



Proof 1. From propositionjSj the solution for small e consists of two rarefaction 
waves and it can be easily checked that the limit intermediate density is zero 
(thanks to proposition [l] point 3, ^/e[P{p)\i = Ole"^)). Since pi tends to p*, 
we have limAl = — oo and limAl — U£. The limit of the 2-rarefaction wave is 
a contact wave (since limA'!i_ ^ limA"^ = Uj.). Let us look at the limit of the 
1-rarefaction wave. 

For each possible speed s G [Al , ug] of the rarefaction wave connecting the 
left state to the intermediate state, we have: 

s = X.{p{s),q{s)) = ^ - Vep'ipis)). (68) 

The state {p{s),q{s)) belongs to the integral curve and then for p{s) < pi 
and according to proposition [l] point 3, q{s)/p{s) tends to ug. If s 7^ w^, 



equation (68 1 implies that p{s) has to tend to p* and that \im ep{p{s)) — since 
\im ep'{p{s)) is finite. This yields the definition of a declustering wave, given in 
definition [T] 

2. We have lim/il(p,.) = prUg > Qr- According to proposition [sj in order to 
have a solution which is the limit of two shock waves, we have to prescribe: 

h%{pi) = PiUr + V£\ —\/[p]r\p{p)Yr < 9£ = PiUt, 
\ Pr 

and in the limit: 




/pi < Ui. 

PrP* 

If h'^^{p() > Qi, the solution is the limit of a 1-rarefaction wave and a 2-shock 
wave. According to this discussion, we have to consider two cases but we will 
see that the limits of the two cases are the same. 

The first case corresponds to the limit of two shock waves. The intermediate 
density is then greater than the left one and so tends to p* too. Besides we 
have: 



pue ' 

and then by dividing by p, we get 



-\/[p]e[^Pip)]e = P'^r + \\ —\/[p\r[ep{p)]r, (69) 

Pi V Pr 



[PWP{P)V , [P\r[£p{p)]r ^ 

V PPe V PP"- 

The last inequality implies that the limit of sp{p) is finite. We denote it by p. 



Taking the limit in the equality ( 69 ) , we obtain 




Ui^Ur + \ v^. 

P Pr 
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Hence we find the value of p as in proposition [5] The propagation speeds are 



easily deduced from the limit e — 7> in ( 67 1 



If the solution is the limit of a combination of a rarefaction wave and a shock 
wave, then the intermediate states satisfies. 



pug - py^[P{p)]l = pUr + A —^/[p]r[ep{p)]r. 

y Pr 

If the intermediate density did tend to p* , then we would obtain ui — Ur 
(thanks to proposition [T]) , which is impossible. Thus the intermediate density 
tends to p* and the previous expression tends to: 



p*Ue = p*Ur + \ —\/{p* - Pr)P, 
V Pr 



which yields the expected result. Since the pressure is positive, Al tends to — oo 
which implies that the rarefaction wave turns into a shock wave with infinite 
propagation speed, i.e. a declustering wave. 

3. From proposition [3] the solution is the limit of a combination of a 1- 
rarefaction wave and a 2-shock wave. Using the fact that pui = pUr-, the inter- 
mediate state (p, q) satisfies: 

ps[P{p)]e = .[^^[pUep{p)]r. 

V Pi 

So [p]r[ep{p)]r (since eP{p) < eP{pi) tends to zero) and either p tends to 
Pr or sp{p) tends to zero. Actually, whatever the limit value of the intermedi- 
ate density, the intermediate state disappears. Indeed, let us consider all the 
possible cases. 

Either p p^. Then the 2-shock wave disappears and the 1-rarefaction 
wave tends to the sum of a declustering wave and a contact wave, which can be 
proven as in the case 1 of this proof. 

Or limp £]pr,p*[. It is easy to check that the 2-shock wave becomes a 
contact wave and the intermediate state disappears since A^, Al — >■ Uf. Like in 
the case 1 of this proof, the 1-rarefaction wave leads to declustering wave and a 
contact wave, which superimposes on the one coming from the 2-shock wave. 

Or lim p = p* . Then we know that ep{p) tends to zero. Thus the 1- 
rarefaction wave tends to a declustering wave and the 2-shock wave tends to 
a contact wave. 

In all the cases, the limit solution is the one given in proposition [5] 

Finally, we consider the case where both the left and right asymptotic states 
are congested: pi ~ pr = p* ■ Besides, we assume that pg and Pr are finite. 

Proposition 6 (case pg = p* , p^ = p* , Pi > Pr) 

1. Ifui^Ur < 0, then the solution consists of one declustering wave connect- 
ing the left state {p*,qe,pe) to a congested and pressureless state {p*,qi,0), 
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then two contact waves connecting {p*,qi,0) to vacuum and then vacuum 
to {p*,qr,0) and another declustering wave connecting {p*,qr,0) to the 
right state {pr,qr,0): 

, * _ s declust. , i „, contact „, vacuum ,„ contact , » „, declust. , » _ , 

(P ,qe,Pe) — 5- (p ,qe,0) — > (0,qe,0) — > (0, gr,0) — > (p ,(jr,0) — > (p ,<jr,Pr). 

2. If U£ — Ur > 0, then the solution consists of two shock waves with infinite 
propagation speed connecting the left state {p*,q£,p£) to an intermediate 
congested state {p* , q, +oo) with infinite pressure and then this intermedi- 
ate state to the right state {p*,qr,Pr)-' 

/ _ \ shock / - , \ shock / * — \ 

(P ,qe,Pe) — > [P ,q,+oo) — > {p ,qr,Prj, 

where q is the unique solution of: 

[qli ^ fPr 
[q]r \Pe 

3. Iful — Ur, then the solution consists of a uniform constant state {p* , qr,pr)- 

Proof 1. The arguments are similar to those used in the proof of the first case 
of the previous proposition. 

2. From proposition [3] the solution is the limit of two shock waves if for 
small e, we have the following inequalities: 



KiPi) = Pt^r + ^/£\j'^\J[pYr[p{p)Yr < = Pt^l, 
y Pr 

ht{Pr) = Plui + VeJ ^\/[pYr\pip)Vr > = pt-Ur- 

y Pe 

Here, these inequalities are always satisfied : their limit is Ur < ui since 
is bounded and [pjf 0. The intermediate density p satisfies: 



[pWp{p)V , [P]r[ep{p)]r 
Ui — \ ; — Ur + \ ; . 

V PPe V PP^ 

Therefore, sp{p) cannot be bounded (which would imply ug = Ur since [p]i and 
[p]r tends to zero). So ep{p) tends to +00 as e 0. From the Rankine-Hugoniot 
relations, we have 

[pq + ep{p)]e[p]e = [q]'j, 
[pq + ep{p)]r[p]r = [q]l 
Taking the limit of their quotient, we obtain: 

limM*^=(Mi)', (70) 

[£P{p)]r[p]r \[q]rj 
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Besides, we have: 

[ep{p)]r[p]r -^^0 [p]r p* ~ p^ 

where the last equivalence results from the fact that {p* — p) — o{e^) and 
{p* — pi^r) — 0{e^). Finally, we have: 

P* - Pi _ ( epip, 



P* - Pr \£piPe) 



This last result, combined with eq. (70), provides an equation for the interme- 
diate momentum. 

3. Let us suppose that pf > p^. The intermediate density satisfies pf. < 
yO^ < pI and so tends to p* . The intermediate momentum q tends to p*Ui. The 
1-rarefaction wave tends to a shock wave with infinite propagation speed (since 
Al and Al tend to — oo), i.e. a declustering wave. We have now to determine 
p. We have 

J-V[p]r[ep{p)]r = -V's[P{p)]e. 

V Pr 

Since p- -> 0, we have -[P{p)]i ^ {pi - p)P'{pi). We have P'{pt) ^ C(p* - 
Pe)~'^ . Then we have P'{pi) ^ C(e~)^T^ = Ce^^^^ (since (p* — pg) = 
0{e^)). Besides {pi — p) < (p* — Pr) — 0{e^) and [p\r — 0{e^) and , so we 
have 

P [P\r 

So sp{p) tends to Pr- The 2-shock wave disappears (since the intermediate and 
the right state are identical). Then the limit solution consists of an instanta- 
neous propagation of the right state. 



A. 2 The one-dimensional cluster collisions 

The Riemann problem where both the left and right states are congested is 
ill-posed since infinite pressure may appear to correct the discontinuity in the 
pressure in the incompressible domain. Like in fl7 |ll2l [5]. we have to restrict to 
the collision of finite congested domains. Consider two one dimensional clusters 
which collides at a time t^- Before collision, the left (resp. right) cluster at time 
t < t(. extends between ai{t) and bi{t) (resp. aj.{t) and br(t)) and moves with 
speed: 

U( — a'i{t) = b'^{t) (resp. Ur — a'j.{t) — b'j.{t)). 

After collision, the two clusters aggregate and form a new cluster at time t > tc 
extending between a{t) and b{t) and moving with speed u — a'{t) — b'(t). 
Therefore, p and u are given for t < tc by: 

P = P*'^[ai{t),beit)] + P*Ma^(t),br.it)], U = Ull[ai,(t),bi.{t)\ + Wrl[a,(t),6,(t)] , 
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and for t > tc by: 

P = P*l[a(t),6(t)], = ul[a(t)Xt)]- 

where 1/ denotes the indicator function of interval /. We denote by to = 
bi{tc) — ar(tc) the cohision point. We look for a pressure written as p{x,t) = 
Tr{x)d{t — tc). Conditions to have such kind of solution for the one-dimensional 
version of (I8|-(l9])-(ll0| was obtained in [1^. We report them in the following 
proposition. 

Proposition 7 1- Supposing thatp{x,t) ~ TT{x)S{t — tc) where tt is continuous 
and zero outside the clusters, then u and tt have to satisfy 



{u — ui){m — a{tc)) + (m — Ur){h{tc) — to) = 0, 



■k{x) 



p* {u — ui){m — x) 

+p*{u - Ur){b{tc) - m), 

p*{u - Ur){b{tc) - X), 



ifxe [a{tc),m] 
ifx£ [m,b{tc)]: 



(71) 

(72) 



2 - Under conditions ( 71)-{ 72), {p,pu,p) is a solution (in a distributional sense) 
to the one- dimensional version of system 



10). 



Remark : We note that the solutions where the two clusters aggregate after the 
collision is one among possibly multiple solutions. Another solution is given by 
a bounce of the two clusters one against each other leading to a reflection of 
the velocities. There are infinitely many such reflections which preserve total 
momentum. 



B Appendix : The two dimensional full time 
and space discretization 

The Direct method: We consider the 2D case with domain Q, — [0, 1] x 
[0, 1]. Denote {x^, y.j) = {iAxjAy), z = 0, • • • , Mi; j = 0, • • • , M2, where Mi = 
l/Ax,M2 = 1/Ay. Let U = {p,q)^,q = (91,92)^ and f/,j = U(xj,yj). To 
simplify the notation, wc define 

Then ^ and ^ can be written as 

dtp + d^qi + dyq2 = 0, 

dtq + d^F{U) + dyG{U) + V,(epi(p)) = 0. 

Denote 

pn_('-^+epoipn\ ^._( ] 
[ ^# j' 'W+epoipnJ' 
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%. 

u defined as follows 



where D'f ,jU, D\ -u are the centered difference operators for any scalar functions 



u • 



We also define the eigenvalues of the Jacobian matrix for two one-dimensional 
hyperbolic system as follows: 



P P ^ P P ^ 

With the above notations, the full discretization of the scheme takes the follow- 
ing form: 





■ 






At 
















At 





J / Ay \ ~^'3- 
I ''Ayi^'i'^+i 



where the fluxes are 



and 



= max{|Ag|, |A« J, |Ag|, |a|^\J}, 
=max{|AgUAg+J,|AgUA|5+J}. 

Similar to the one-dimensional case, by inserting the momentum equation 
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into the density equation, we can get the following elliptic equation 



„n+l 



At^ r 1 



4 Aa;2 



At^ ( 1 

1 



+ 



AxAy 
1 



Ay 



At 
2Ax 

At 
2Ay 



iF7+i/2,j+i f^ - iF'U/2,j+i f^ - iF7+l/2,^l f^ + {FU/2,j-l)^'^ 

(73) 



Here (-^"+1/2 j)*^^' is the first component of the vector .^"+1/2 j- Also similar to 
the one-dimensional case, we solve the above elliptic equation to get first Pi~^^ 
then And once p"+^ is derived, we can get 5"+^ explicitly. 



n+l 



At 



At 
Ay 



G 



G 



AtV,j(epi) 



n+l 



'^■J Ax ' 

Gauge method: The Gauge method can be implemented in a similar way 
as the ID case. Indeed, we have the same elliptic equation for pi (731 and the 
following equations: 

1 r n+l r, n+l I n+l 1 , If n+l r, n+l , n+l 1 



4Aa;2 



At 



2Ax 



^i+h-i(Pi+i,i ~ Pi,i) ~ ^i-l,iiPi,i ~ P'i-l.j) 



2Ay 



^i,j+\{Pi,j + l ~ Pi.j) ~ (^i,j-^{Pi,j ~ Pi,j~-l) 



4Aa;2 



•■■J 

-I — 
2\ Aa;2 



4Ay2 »j+2 2fij 



■n+l I j 



(f: 



i+3/2j 



)''^ - iF7+i/2,,Y'^ - {FU/2j'^ + {F-_^/2,r^ 
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AxAy 
1 



1 1 / \ 1 / \ 



This is the Gauge 2 method. And similar to the one dimensional case, we will 
mainly test the Gauge 1 method with a smaller stencil in the Laplace equat 
of f and P: 



;ion 



4- 



J-lJ 



At 



2Ax 



^ fpn+l 9 pn+1 I pn+1 ] , ^ [pn+l _ opn+l i pn+1 1 



+ 



1 

AxAy 
1 



AxAy 
1 

A? 



(-^"+l/2j + l)''^' - (-P'"-l/2j + l)''^' - (-PTh 



+i/2.-i)(^' + (i^r-i/2,-i)^^) 

Different from the one-dimensional 



case, now we are facing three elliptic 

equations in each time step. In section [5. 2[ we will see some simulation results 
and comparison for these two methods. 
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